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Abstract 

Tidal Disruption of stars by supermassive central black holes from dense rotating star clusters is 
modelled by high-accuracy direct Af-body simulation. As in a previous paper on spherical star clusters 
we study the time evolution of the stellar tidal disruption rate and the origin of tidally disrupted stars, 
now according to several classes of orbits which only occur in axisymmetric systems (short axis tube 
and saucer). Compared with that in spherical systems, we found a higher TD rate in axisymmetric 
systems. The enhancement can be explained by an enlarged loss-cone in phase space which is raised 
from the fact that total angular momentum J is not conserved. As in the case of spherical systems, the 
distribution of the last apocenter distance of tidally accreted stars peaks at the classical critical radius. 
However, the angular distribution of the origin of the accreted stars reveals interesting features. Inside 
the influence radius of the supermassive black hole the angular distribution of disrupted stars has a 
conspicuous bimodal structure with a local minimum near the equatorial plane. Outside the influence 
radius this dependence is weak. We show that the bimodal structure of orbital parameters can be 
explained by the presence of two families of regular orbits, namely short axis tube and saucer orbits. 

Also the consequences of our results for the loss cone in axisymmetric galactic nuclei are presented. 

Subject headings: black holes - galactic nuclei - stellar dynamics 


1. INTRODUCTION 

A large fraction of galaxies show evidence of supermas¬ 
sive black holes (henceforth SMBH) residing in their cen¬ 
ter. They are typically embedded in nuclear star clusters 
(NSC); if resolution allows to observe the NSCs, they are 
among the densest clusters known. Their size is similar 
to galactic globul ar clusters, but t h ey ar e much heav¬ 
ier and brighter (|Boker et al.l I2002L 120041) . In massive 
galaxies NSCs may not be significant or even do not ex¬ 
ist, however, the SMBHs are still surrounded by enor¬ 
mous number of stars. SMBH residing in these NSCs 
will tidally disrupt stars that come close to its tidal ra¬ 
dius and eventually accrete the gaseous debris, whic h can 
light up the central SMBH fo r a period of time (|Reesl 
119881 lEvans fc Kochan^ll989[) . This kind of event is a 
useful tool to examine the relativistic physics near SMBH 
since the disruption occurs at a place very close to the 
BH’s Schwarzschild radius. Also it can help us to in¬ 
vestigate SMBH in non-active galactic center. Although 
tidal disruption of stars has been proposed for almost 
half a century, only until last decade do people realize 
the importance of such events, after th e discovery of a 
dozens of tidal disruptio n candidate s (iKomossal 1200^ 
iKomossa fc Merritt 1200811 . iLiu et all (|2014D discovered 
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a candidate of binary SMBH system by analyzing the 
break in the light curve of TD event, demonstrate it 
as a promising tool for searching hidden SMBH bina¬ 
ries in quiescent galactic center. In order to compute the 
tidal disruption event rate, many t heoretic works have 
been done in the past few decades (iFrank fc Rees 1976; 
Lightman fc ShaDirolll977l: iMagorrian fc TremaindllQ^ : 
Wang fc Merrittll2Q04[) . The core of the story is loss cone 
theory, which was first established in the case of spherical 
symmetric systems. 

Stars with orbital pericenter smaller than the tidal ra¬ 
dius vt are defined to be inside the loss cone, with rt be 
expressed by 
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n = OLV^ - - 

b — n 
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where r*, m* are the radius and mass of a star, n 
is its polytropic index (assuming the stellar structure 
can be approximated by a polytropic sphere) and a is 
a free parameter used by us for scaling. Stars with an¬ 
gular momentum J < Ji^ ylGM^rt are inside the 
loss cone. Typically, loss cone stars are consumed in 
dynamical time scales. If no new star is supplied to loss 
cone, there will be no more tidal disruption event. Based 
on the status of loss cone, it can be divided into two 
regime, namely empty and full loss cone. Due to the 
short “lifetime” of the loss cone stars, loss cone will be¬ 
come empty quickly. Refilling of loss cone happens in 
relaxation timescales and is often referred to as diffusion 
process in angular momentum space. Thus in empty loss 
cone regime it is the refilling rate which controls the dis¬ 
ruption rate. Note that throughout this paper, and like 
in most if not all of the cited papers on tidal accretion 
of stars onto SMBH, we assume that a star is disrupted 
completely at rt and all its mass, energy and angular 
momentum absorbed momentarily by the SMBH. We 
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know that this is not realistic, and more detailed nu¬ 
merical models of the process of disruption, possible disk 
formation and accretion show that only fractions of the 
materi al are absorbed into the SMBH after a number of 
orbit s (jGuillochon fc Ramirez-Rnidl2013HHavasaki et alJ 
1201,111 . However, the assumption that the process is fast 
is reasonable compared to the orbital time scales of stars 
further out in the cluster^ _ 

In a previous paper (|Zhong et al.l (|2014ll . henceforth 
Paper I) we have shown that the classical loss cone ap¬ 
proximation, for a spherically symmetric system in the 
diffusive empty loss cone regime, can be well reproduced 
by large direct iV-body models with tidal accretion of 
stars onto SMBH. Now we are focusing on the general¬ 
ization to axisymmetric galactic nuclei and compare our 
new results in an otherwise very similar study to those 
of Paper I. 

Tidal disruption of stars is one possible way of growth 
for SMBH, especially in quiescent galactic nuclei. Since 
most models assumed spherical stellar clusters, SMBH 
growth rates by tidal disruption are very low, limited by 
the very long relaxation time to refill the loss cone, and 
the contribution of the process to the overall growth of 
SMBH is considered as relatively insignificant. However, 
the stellar distribution in real galactic nuclei might not 
be spherically symmetric. Many galactic nuclei show ev¬ 
idence o f rotation in their centers, even very close to the 
SMBH (IMivoshi et al.l 119951 INeufeld & MalonevI 119951 
Inreenhil1eral.lll99lb ^ 

According to the current standard model of structure 
formation massive galaxies have undergone quite signif¬ 
icant mergers (in number and mass ratio). Numerical 
models of the merging process of galaxies show that 
the merger remnant shows rotation, axial symmetry or 
even triaxiality in the central regions dlHian et al1l2011l: 
Preto et al.ll201 ll iGualandris fc MerrittJl2012tlBois et al l 

2nTl|b 

In the center of our own Milky Way the NSC can be 
observed in unparalleled high resolution (|Feldmeier et al.l 
12014 iSchodel et al.l[20l^ . It consists of 1.4 x IO^Mq 
within its effective radius (4.2 pc); kin ematic data in¬ 
dicat e that it possesses bulk rotation (|Feldmeier et al.l 
12014(1 . The formation mechanism of NSCs is still un¬ 
der debate. There are two scenarios, in situ forma¬ 
tion (|Milosavlievicl[2?)?)3 l and a sinking sc enario (globular 
clust e r sink to the cen ter and merge) (I Tremaine et al.l 
Il975t iLotz et all l2001h. NSCs in a sample of nearby 
galaxies observed bv iSeth et al.l (|2006L 1200811 show that 
these objects are non-spherical and even contain multi- 
component (younger disk plus older spherical com- 
ponent),which favor the in situ scenario. However, 
lAntonini et al.l (|2012ll have performed a series of A^-body 
simulations to study the formation of NSCs, which sup¬ 
port the sinking scenario. The model NSCs formed in 
their simulations by merging between infalling globular 
clusters initially have mildly triaxial shape. After the 
final infall, the shape of the NSC will gradually become 
axisymmetric in following dynamical evolution. 

Despite the debate between different formation scenar¬ 
ios, we think that it is quite likely that NSC are non- 
spherical. This provides a good motivation to study the 
tidal disruption rate in axisymmetric (and triaxial) clus¬ 
ters. Some works have already been done but the mission 


is not over. iFiestas fc Spurzeinl (2010 1 used 2D Fokker - 
Plank model ( Einsel fc Spurzeia 199S ; iKim et ^ l2002f) 
to study rotating dense stellar clu sters with BHs an d 
cross checked with iV-body models ([Fiestas et al.ll201^ . 
Both works find that BH embedded in rotating model 
have higher tidal disruption rate (hereafter TDR) com¬ 
pare to spherical models. BH mass at the end of simu¬ 
lation is roughly 20% higher in rotating case. They hnd 
an excess of accreted prograde rotating stars which are 
originated mainly outside the influence radius Vh and call 
for a further investigation of the roles of stars with non- 
conserve d Jx,Jv angular momen tum. As shown by the 
works of iMerritt &: Poonl 1)200411 , in non-spherical sys¬ 
tems chaotic orbits (existing in regions outside Vh) can 
keep the loss cone full for sufficient long time, thus tidal 
disruption can contribute a lot of mass within Hubble 
time and could play an important role in the BH growth 
across cosmic time. 

On the other hand, t he loss cone itself might be en - 
larged as pointed out bv lMaeorrian fc Tremaind (I1999I1 . 
due to the fact that angular mome ntum J is not con¬ 
served in axisymmetric potential. IVasiliev fc Merrittl 
(|2013ll confirmed this picture in a detailed analysis of the 
loss cone problem in axisymmetric galactic nuclei. They 
analyzed the depletion and refilling of loss cone orbits 
and found that tidal disruption rates could be increased 
by a moderate factor due to axisymmetry as compared 
to spherical symmetry. In their work chaotic orbits with 
low angular momentum, which can reach just outside 
the influence radius at apocenter, but also get close to 
the central SMBH at pericenter, cause some difficulty in 
comparis on with Fokker - Planc k models, as was already 
found by iMalkov et al.l (|1993ll (Note that the last au¬ 
thor of this paper is t he sa me person than the last au¬ 
thor of IMalkov et al.l (119931) . there was a mistake in re¬ 
translating the name from Russian language). 

In this work we follow an experimental numerical ap¬ 
proach to the problem, following Paper I for the case of 
spherically symmetric systems. We treat particle num¬ 
ber and tidal radius as free parameters and analyze the 
tidal accretion rate of the system as a function of the 
strength of deviation from spherical symmetry. We mea¬ 
sure the shape of the loss cone in axisymmetric potential 
and and characterize the characteristic orbits of stars in 
the loss cone. We find that it is indeed enlarged and can 
account for the higher TDR as compared to spherically 
symmetric galactic nuclei. 

This paper is organized as follows: we describe the 
model setup of the simulation in Section 2 and present 
the result of TDR measurement in Section 3. Section 
4 is devoted to the measurement of loss cone shape in 
axisymmetric potential and we demonstrate the enlarge¬ 
ment of loss cone. In Section 5, we present the result for 
the origin and orbital classification of disrupted stars. 
In Section 6, we discuss the potential application of our 
results. 


2. WBODY MODEL 

We adopt the st a ndard A^-body unit definitions from 
iHeggie &: Mathi^ (|1986D . namely G = M = 1 and E 
= —1/4, where G is the gravitational constant, M is the 
total mass of the model cluster and E is the total energy. 
In our A^-body models we assume that all the particles 
have the same mass, so m = 1 /N, where m is the particle 
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mass and N is the total particle number. To preserve the 
scale invariance of our iV-body simulations we fix the 
initial black hole mass relative to the total mass of the 
star cluster (0.01) and use the particle number and the 
tidal radius in A^-body units (which is a dimensionless 
number) as free parameters. We have shown in Paper I 
that the method of scaling to realistic parameters for N 
and rt can be used to obtain astrophysically meaningful 
results from the collection of our models. In order to 
support our scaling procedure we even do not change 
the tidal radius during the simulations - since the BH 
mass changes within one order of magnitude only during 
the simulation, relative changes in tidal radius are small 
(notice that r* oc (M,)^/^). 

The initial distribution of particles follows a general- 
ize d King model with rotation. The distributi on function 
is (lEinsel fc Spurzeinlll999l : lErnst et al.ll2007ll 

f(E, J^)=C ■ [exp(—^) - 1] • exp(-^^), (2) 

where gk is the King velocity dispersion and Do is a 
characteristic angular velocity. Since we are considering 
an isolated system, the is set to 0. This rotating King 
model has two dimensionless parameters: Wq and wq- 
The King parameter Wq = where $o is the cen¬ 

tral potential, controls the degree of central concentra¬ 
tion. And the rotation parameter wq = \/9/(47rGpo)’Do, 
where po is the central density, controls the degree of 
rotation, wq = 0 will reduce the model to a usual non¬ 
rotating spherically symmetric King model. 

We limit our current study to only one concentra¬ 
tion parameter Wo = 6 and two rotation parameters 
Wo = 0.3,0.6; the density profile of King model with this 
concentration is similar to that of the Plummer model 
used in Paper I, so it is possible to compare with the 
previous results and focus on the effects of rotation and 
axial symmetry on l y. Th e rotation is moderate (cf. e.g. 
lEinsel fc Smirzenil {199^) and resembles that of Milky 
Way globular clusters. 

For completeness we also employ non-rotating King 
model with Wq = 6 and wq = 0.0, which is used as a 
fiducial model and also a bridge to the results of Paper I, 
conhrming our claim that it indeed closely resembles the 
results for the Plummer model used in Paper I (e.g. in 
the evolution of the TDR). In another test run we used a 
larger rotation with wq = 0.9 - it experienced an unstable 
stage during which a bar formed but quickly disappeared. 
This bar formation could probably be identif i ed wit h the 
radial orbit instability of I Aguilar fc Merritt! (|I990D . We 
note that our standard models with wq = 0.3, 0.6 remain 
fully axisymmetric during the entire simulation; to study 
tidal disruption in triaxial systems with bars is beyond 
the scope of our current paper. 

Fig. [T] shows the axial ratio (c/a) of the model clusters 
as a function of radius up to r = 2.0 (within which most 
of stars are located). We estimate the axial ratio for 
both rotating models, using the moment of inertia tensor 
measured in concentric shells. One can see c/a is close to 
I at the innermost part and decreases outward: ujq = 0.3 
model decreases slowly to its minimum value 0.9; wq = 
0.6 model decreases faster and has a minimum value 0.71. 
If we measure the c/a for the whole cluster, the results for 


the two models are 0.9 (wq = 0.3) and 0.75 (wq = 0.6). 
Fig. [T] also shows that c/a is almost unchanged during 
long time evolution, except for the inner part of ujq = 0.6 
model, which exhibits slight decrease. 



r [NB] 

Figure 1. Axial ratio for rotating models as a function of radius. 
For each model we show the axial ratio measured at different evo¬ 
lution stage: T=0 (red); T=500 (green); T=1000 (blue). Lines 
with symbols are results for ojq = 0.3 model, lines without symbols 
are for ljq = 0.6 models. 

In rotating systems, there is a phenomenon called 
gravo-gyro instability, which is caused by the nega¬ 
tive specihc moment of inertia (jinagaki fc Hachisnll19781 : 
iHachisul 11979111982D . This kind of instability happens 
in long term evolution of rotating cluster which is m uch 
longer than our integration time (lErnst et al.ll2007t) . 

The model set is summarized in Table [T] 


Table 1 

Full set of our model runs. 


Model 

N/K 

iOO 

rt 

T 

R20w00 

64 

0.0 

10“^ 

1500 

R30w00 

128 

0.0 

10-3 

1600 

R21w00 

64 

0.0 

lO-'^ 

1500 

R31w00 

128 

0.0 

10-4 

1300 

R20w03 

64 

0.3 

1 

O 

1500 

R30w03 

128 

0.3 

10-3 

1500 

R21w03 

64 

0.3 

10-4 

2600 

R31w03 

128 

0.3 

10-4 

2000 

R20w06 

64 

0.6 

1 

O 

1500 

R30w06 

128 

0.6 

10-3 

1500 

R21w06 

64 

0.6 

10-4 

1600 

R31w06 

128 

0.6 

10-4 

2000 


Note. — Column 1 : Model codename. Column 2 : Particle 
number in the unit of K(=1024). Column 3 : dimensionless rota¬ 
tion parameter. Column 4 : black hole’s tidal radius. Column 5 : 
total integration time, rt and T are in model unit. 


We run the simulation for more than one initial half¬ 
mass relaxation time (fcii); which is estimated using the 
same formula in Paper I and the values can be found 
there as well (Table 2). 

All simulations ar e running with the t/sGRAPE code 
(|Berczik et al.l[2MTI) . which runs with high performance 
(up to 350 Gflop/s per GPU) on our GPU clusters in 
Ileijing (NAOG/CAS). This code is a direct A^-body sim- 
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ulation package, with a high order Hermite integration 
scheme and individual block time steps. A direct N- 
body code evaluates in principle all pairwise forces be¬ 
tween the gravitating particles, and its computational 
complexity per crossing time scales asymptotically with 
iV^; however, it is not to be confused with a simple brute 
force shared time step code, due to the block time steps. 
We refer more interested readers to a general discus- 
sion about iV-body co des and their implementation in 
ISpurzem et al.l (|2011al lbll . The present code is well tested 
and already used to obtain important result s in our ear¬ 
lier la rge scale few million body simulation (|Khan et al.l 

[Ml). 

3. TIDAL DISRUPTION RATE (TDR) 

3.1. Results of our work 

In this section, we present the TDR measured in sim¬ 
ulations with our rotating King models and compare it 
with the TDR of the non-rotating model of Paper I. In 
Fig. 121 we show the TDR (both in terms of mass and 
particle number) as it evolves with time for two differ¬ 
ent tidal radii; in each panel two different rotation pa¬ 
rameters are plotted together with the data of the non¬ 
rotating system. The time is given in units of initial half 
mass relaxation time trh, which is convenient for com¬ 
parison of simulations with different particle numbers. 
To smooth out fluctuations due to particle noise we have 
plotted in the figure the TDR averaged over a time in¬ 
terval (here l/4U?i)- 


r^.=10 ^ r^=10 ^ 



0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 



0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 


t / t / 

Figure 2. TDR as a function of time in units of initial half mass 
relaxation time; averaged over intervals of l/4£^^. Top panels: 
mass accretion rate; bottom panels: particle accretion rate; left 
and right panels for two different tidal radii as indicated. Curves 
with symbols are stand for 128K models, those without symbols 
are stand for 64K models. Line thickness indicate different rotating 
parameters. 

The TDR with a large tidal radius (i.e. rt = 10“^) 
initially quickly rises in the N = QAK model to its peak 
value, and then decreases; for the N = I28Ar model the 
TDR almost decreases from the beginning. The initial 
phase is connected with the formation of a central density 
cusp in the surrounding stellar system and with the pro¬ 
cess of transition from initially full to empty loss cone. 
The BH gains mass from the accreted stars, thus the 
mass ratio between stars and the BH (7 := m/M,) de¬ 


creases with time, and as a result the BH’s random mo¬ 
tion damps. We have discussed in Paper I that the sta¬ 
tus of the loss-cone is connected with the BH’s Brownian 
motion in the sense that once the amplitude of Brownian 
motion is smaller than lOr^ the system enters the empty 
loss-cone regime, during which the cusp and central den¬ 
sity are still growing but TDR begins to fall. In the 
N = 128Ar model, the mass ratio 7 is smaller, so the ini¬ 
tial loss cone depletion is very short, practically invisible 
in the plots, and the subsequent evolution is determined 
by cusp formation and damping of BH motion. 

In the models with small tidal radius (rt = 10“"^), there 
is always an initial growth phase of TDR, followed by the 
convergent approach to a stationary state. Due to the 
small rt their BH growth is slow, thus they need more 
time to achieve the mass required to limit their Brownian 
motion. 

Fig. 12] also shows the TDR dependence on rotation pa¬ 
rameter ujQ as a new result compared to Paper I. For large 
tidal radius {rt = 10 “^), faster rotation will result in a 
higher TDR, note that these models are in empty loss- 
cone regime. Table |2I list out the numbers for TDR mea¬ 
surement. One can see ujq = 0.3 model has a TDR on av¬ 
erage 13 percent higher than wq = 0.0 model. And TDR 
in Wo = 0.6 model is on average 35 percent higher than 
that in wq = 0.0 model. BH mass of these 3 models mea¬ 
sured at T = 1500 are 0.131, 0.143 and 0.167. The frac¬ 
tional increase of final BH mass with incr easing degree 
of rot ation is consistent with the result of [Fiestas et al.l 
(|2012ll . The reason for this dependence of wq is that in 
these systems the effective loss-cone is larger than clas¬ 
sic one in spherical system. We will investigate such 
an enlarged loss-cone in more detail in the next section. 
For small tidal radius {rt = 10“^), however, we observe 
a different behavior of TDR. From beginning to about 
~ 1.5trh, faster rotation result in a smaller TDR! The 
argument presented by iMagorrian fc Tremaine! (|1999f ) 
may provide some hints: if BH’s wandering time-scale 
is shorter than dynamical time-scale, a decrease in TDR 
will happen. We note in the simulation at this early stage 
the BH is quickly wandering due to its small mass and 
slow growth. Furthermore, in axisymmetric systems a 
star’s pericenter distance changes with time (even ignor¬ 
ing irregular perturbations from other stars). So when 
the BH comes back to the place where it was, it may still 
miss the star which is supposed to be disrupted shortly 
before. 


Table 2 

TDR results for rt = 10models. 



No 

No 

No/No 

No 

Ne/No 

0.25 

14.96 

15.41 

1.03 

16.30 

1.09 

0.50 

12.89 

13.67 

1.06 

17.30 

1.34 

0.75 

10.44 

12.07 

1.16 

14.52 

1.39 

1.00 

8.73 

10.32 

1.18 

12.65 

1.45 

1.25 

7.97 

9.73 

1.22 

11.09 

1.39 

1.50 

6.99 

7.93 

1.13 

10.04 

1.44 


Note. — Measured TDR for models with same N (128K) and 
rt (10“®) but different rotating parameters, at different evolution 
time. No is TDR in classic King model (wq = 0); No and Nq are 
results for ojq = 0.3 and ojq = 0.6 models. We also give the boost 
factor No/No and Nq/Nq- 
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Afterwards the system begins to enter the empty loss 
cone regime, and all TDR curves converge to each other; 
for small tidal radius more tidally disrupted stars origi¬ 
nate from inside the BH influence radius, where the sys¬ 
tem is approximately spherically symmetric. Any devia¬ 
tion from spherical symmetry in our rotating models pre¬ 
vails near and outside the influence radius. Convergence 
of TDR reflects the original results obtained in Paper I 
for spherical systems. 



t / trt 


t / trj 


accretion radius vt has been chosen very large compared 
to any realistic tidal radius (for smaller N simulations 
it has to be done in order to get any meaningful results 
on TDR). So, we also have to discuss how to scale down 
the TDR from our simulated values of rt (10“^, 10“^) to 
the small more realistic regime of rt (10“^). This can be 
done by applying scaling relations from known scaling 
laws (obtained e.g. from Paper I and other literature) 
for N, and an empirically determined one for rt. We 
have shown above that the TDR in King (Wg = 6) and 
Plummer models is very similar, so we expect the scaling 
formula (AlO) derived in Paper I for a Plummer model 
to be also valid for our King model used here. So, we 
apply the same boost factor of TDR with respect to the 
axisymmetry of a galactic nucleus for the real galactic 
nucleus as we find here in this paper for our simulated 
systems. For example, in Paper I we estimated the TDR 
of the Milky Way SMBH to be 1.09 x 10“®yr“^ after a 
scaling procedure with respect to N and r*. By fitting 
surface brightness profile to mid-i nfrared images of the 
nuclear cluster in our Milky Way, iSchodel et al.l l)2014[) 
reported the mean ratio between minor and major axes is 
0.71, which is close to our rotating wg = 0.6 model. For 
this model we find a boost of TDR by 35% in our simu¬ 
lations, and we apply the same factor here for the case of 
axisymmetry, to get a higher TDR of 1.47 x 10“^yr“^. 


Figure 3. x axis is time in unit of initial half-mass relaxation 
time, y axis for panel a) and b) is the averaged mass accretion 
rate in given time range (i.e. 1/4 trh)\ V for panel c) and d) 
is the number accretion rate. Panel a) and c) show the result for 
rt = 10“^. Panel b) and d) show the result for rt = 10”'^. 

Fig. compares the TDR of classic King model {W = 
6,a;g = 0.0) with that in Plummer model. In rt = 10“^ 
models, except the initial higher accretion rate in King 
model, the two models have similar TDR in following 
evolution. While in the case of rt = 10“^, King model 
have a higher accretion rate during most of the time, 
but later on they gradually come to the same level as 
the Plummer model. The higher rate in King model 
could be explained by the slightly higher density in the 
core region at beginning. In the following evolution of 
rt = 10“^ models, the two models form cusp similar to 
each other so they have roughly same accretion rate. In 
the case of rt = 10“"^, the initial accretion rate ratio 
^King/Npium IS higher than those in r* = 10“^. BH 
inside King cluster growing faster and also the growth 
of cusp, in the following evolution King model always 
have a higher density in the cusp which in return gives 
a higher accretion rate. Only after the BH gain enough 
mass and become a “static” object, the accretion rate 
slowly reaches a maximum and begins to drop afterward. 

Up to this point, all results were presented in model 
units (V-body units). As in Paper I (see Sect. 5 and 
Appendix therein) we will discuss now any conclusions 
which can be made for the case of real galactic nuclei and 
environments from our results. This will be useful for ob¬ 
servational programmes on TDR. To predict the TDR in 
real galactic nuclei, we use the method of scaling. The 
TDR obtained in our simulations has to be scaled up in 
two ways: first from relatively low (TV ^ 10®) to more 
realistic high particle numbers (TV ^ 10®). Second, our 


3.2. Relation to other current papers in the field 

With regard to the enhancement of TDR in axisym- 
metric systems we have shown th at our results are in 
agreement with iFiestas et ahl (1201211: but recently n umer ¬ 
ical simulations pu blished bv iVasiliev &: Merrittl (|2013ll 
and lVasilT^ (j2014H seem to contradict our findings. They 
claimed that the TDR in axisymmeric nuclei can be a 
few tim es larger than in the spherical case. Also lLi et al.l 
(|2014ll analyzed the distribution of stellar orbits in an ax- 
isymmetric galaxy and found that total number of stars 
that can interact with the central SMBH binary is six 
times larger than in the spherical system. In this sub¬ 
section we will discuss why there is such a discrepancy 
to our results - we find a much smaller enhancement of 
TDR in axisymmetric systems. 

The main difference between the cited papers and our 
work is the initial model. In all of the above mentioned 
papers, a flattened Dehnen model is used (their density 
profile, given the parameters they chose, is identical also 
to the Hernquist model). Their models possess a fixed 
axial ratio (c/a = 0.75) throughout the entire cluster and 
an initial central cusp, while our rotating initial model 
has initially a core density distribution in the center, and 
we have a radial variation of c/a from nearly spherical 
(c/a ~ I) to about c/a ~ 0.7 in the outskirts (see Fig.[T]). 
However, in the radius range where most of the disrupted 
stars originate from (c.f. Fig. El), the system deviates 
significantly from spherical symmetry, thus we can con¬ 
firm that the enhancement of TDR is connected with the 
non-spherical geometry. But in the relevant region of our 
ojg = 0.6 deviation from spherical symmetry is less than 
in the other cited papers, which may be an explanation 
for the weaker effect in our case. 

We also notice that even within the cited other pa¬ 
pers there are some discrepancies in the results even for 
models with the same initial density profile. For exam- 
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pi e, the enhancemen t of th e number of accreted stars 
in iVasiliev fc Merritt! (|2013f) was smaller than 100% (see 
Table 2 in their paper), while iLi et al.l (|2014f) found a 
factor of six. On the oth er hand some of the models in 
IVasiliev fc MerrittI (j201,‘l[l only show mild enhancement 
which is in the same level as ours. Another example 
comes from the debate about the “final parsec problem” 
in SM BH binary evolutio n. Based on their simulation 
results. iKhan et all (I2013I1 claimed that the “final parsec 
probl em” is not a pr o blem in axisymmetric host galaxies, 
while IVasiliev et al.l (|20l4l reached an opposite conclu¬ 
sion according to their simulation. We notice that both 
of these work employs similar flattened galaxies model, 
however, they used a different method to generate the 
i nitial model. _ 

IVasiliev fc MerrittI ()201,1f l , IVasilievI (120141 1 and 
IVasiliev et al.l (1201411 utilized the orbital superposition 
method of (lSchwarzschildlll979ll to c onstru ct the ir model. 


On the other hand, iKhan et al.l (|2013ll and ILi et al.l 
(|2014ll used another method called “adiabatic squeeze 
techni que” developed by iHollcv-Bockelmann et al.l 
(|2001tl . We notice that in the process of adiabatic 
squeeze, which contains a step which applies a slow and 
smooth velocity change on the stars in the 2 ; direction. 
This step may artificially reduce the energy and angular 
momentum of the stars in the model cluster. Although 
the radius and velocity vectors of the stars are rescaled 
after the squeeze, it is not clear how the rescaling affects 
the phase space distribution. Thus it might be possible 
that the process produces more stars of low energy 
and low angular momentum. Another evidenc e of a 
similar effect can be derived from IVasilievI (|2014f) : while 
they still use the orbital superposition method they 
changed the generation of their initial model so that it 
creates more low energy and low angular momentum 
stars. In their test run (Fig. 2 in their paper) we see 
a much larger en hancement of the number of accreted 
stars compared to IVasiliev fc Merritd (I2013D . So to add 
more low energy and angular momentum stars seems 
to be promising in abridging the d if ferent enha ncement 
factor s between IVasiliev fc Merritt] (1201,111 and ILi et al.l 
()2014ll . We suggest that a detailed comparison between 
models constructed with these two methods (and their 
phase space distribution) should be performed in order 
to explain the discrepanc y. _ 

According to ILi et al.l (j2014[ l the central two parsecs 
of their model galaxy exhibit a slight triaxiality, which 
could also introduce some additional centrophilic orbits, 
thus increase the number of stars that can interact with 
the central SMBH binary. 

Before finishing this se ction, we want to make a final 
remark on the result of iLi et al.l 1)20141 1 . Their model 
integrates individual orbits in a fixed model potential 
with one SMBH in the center, in a static way. So the 
number of stars that can interact with the central SMBH 
binary according to their results should be considered 
as an upper limit. Once two-body relaxation is turned 
on, some of the stars that are supposed to be inside the 
loss cone might be scattered out. And the presence of a 
SMBH binary in an evolving system may also affect the 
result of how many stars can interact with them. 


4. LOSS CONE IN AXISYMMETRIC POTENTIAL 


Eirst, we summarize the loss cone theory for stellar or¬ 
bits in a spherically symmetric gravitational potential, 
in order to discuss different behavior in an axisymmetric 
potential later. If a stellar orbit has a pericenter dis¬ 
tance less than the tidal radius it is considered to be in 
the loss cone. In spherical symmetry the boundary of 
the loss cone can be expressed in t erms of a critical loss 
cone an gular momentum Jir \/2 GM,rt (if r ^ rt, 
cf. e.g. [Amaro-Seoane et al.l (|2004lH . The loss cone is 
then defined as the region in phase space where the an¬ 
gular momentum J of a star fulfils J < Jic- All stars 
inside the loss cone will reach the tidal radius within 
a dynamical (orbit) time scale. As a consequence the 
loss cone would become empty in that relatively short 
time. Once a star is inside the loss cone and reaches 
the tidal radius, we assume that it will be destroyed by 
the BH’s tidal force instantaneously and add its total 
mass to the black hole at the same moment. Most au¬ 
thors studying stellar dynamics and TDR of star clusters 
around a BH used similar approximations. Rees (1988) 
already argued that the stellar debris after tidal disrup¬ 
tion will make several orbits until it is finally accreted 
by the BH; nevertheless the orbital time near the BH is 
very short compared to the original orbital time of the 
star before its disrup tion. Recent detailed simulations 
on tidal disruptions (iGuillochon fc Ramirez-Ruid 120131 : 
iHavasaki et al.l[2013L i2015|l show that in some case not 
all material of the star may be accreted and that general 
assumptions about the tidal fallback rate are not correct; 
for example in a longer lived accretion disk may form, 
which would delay the black hole growth. In a spherical 
system, without interactions between the stars, angular 
momentum J would be strictly conserved. So, without 
any repopulation of the loss cone, the accretion process 
would stop after a few dynamical times. But stars do in¬ 
teract with each other while moving inside the star clus¬ 
ter by two body relaxation through mutual encounters; 
in this process they can exchange angular momentum 
and energy and so the loss cone will be repopulated in 
the two body relaxation time scale, which is generally 
long c ompared to the dynamical time (jCohn fc KulsrudI 
II978I : lAmaro-Seoane et*^ 12004( 1. The repopulation of 
the loss cone is modelled in these papers as a diffusive 
process using the Eokker-Planck approximation. 

In an axisymmetric potential, the situation is more 
complex since J is not a conserved quantity. It changes 
continuously due to the non-central force resulting from 
the geometry of the potential. In this case, stars with 
J > Jic may have a chance to drift into the loss-cone and 
get disrupted. In other words, the loss cone is enlarged 
in the J dimension in axisymmetric potential. However, 
the 2 component of angular momentum is still con- 
served, so a solid boundary o f the loss-cone is Jz < Jic- 
iMagorrian fc Tremainel (119991 1 inves tigated this topic us¬ 
ing a symplectic map introduced bv lTouma fc: Tremainel 
(|1997ll . In this work, we analyze the enlarged loss cone 
in phase space in terms of energy A, modulus of angular 
momentum J and the 2 ; component of angular momen¬ 
tum Jz for s tellar orbits near the BH . We u se a different 
approach as IMagorrian fc Tremainel (|1999l l here, which 
is based on a numerical particle scattering experiment. 
In what follows, we hrst describe the method we used in 
this experiment, then present our results. 

First step, we need to know the smooth gravi- 
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tational potential as a function of position without 
the fluctuations due to the discrete particle structure. 
We use a so-called self-co nsistent field code (SCF, 
iHernonist fc OstrikeH 1)19921 11 to generate the analyti¬ 
cal function for the gravitational potential. The expan¬ 
sion coefficients Cim , Dim , Eim , Fim. used in computin g 
forces (Eq.(3.21)-(3.23) in iHernauist fc Ostrikerl (|1992ll l 
are computed based on snapshot data generated during 
the direct A^-body simulation. By default, the code uses 
radial basis functions labeled from n = 0 to Umax = 14, 
and spherical harmonic function truncated at Imax = 10 . 

Our particle distribution is self consistently achieved 
as a consequence of the co-evolution of stars and BH. 
Using the SCF code means that all two-body interac¬ 
tions are smoothed out in the experiment. Because 
we assume that most of the two-body interactions hap- 
pe ns during the apocenter p assage, which is also used 
in iTouma &: Tremainel (|1997ll . After getting the coeffi¬ 
cients, we can calculate the acceleration, jerk and do or¬ 
bit integration using a Hermite 4*^ integrator with vari¬ 
able time steps, developed by ourselves. This code works 
very well and the energy and angular momentum errors 
of the test particle stays in the level of 10 “® over long 
time integration. In an axisymmetric system all coeffi¬ 
cient with m 0 should be 0. But in practice one will 
get some small numbers very close to 0 due to particle 
noise. We just ignore these terms, otherwise Jz would 
no longer be conserved. We also ignore coefficients with 
odd I, because the rotating system should be symmet¬ 
ric about the equatorial plane and do not have pear-like 
shape. 

Next step is to generate initial positions and velocities 
for test particles. The basic idea of this experiment is 
to do parameter space scanning. We uniformly sample 
E,J and Jz, all test particles are initially put at their 
apocenter. Firstly, we choose a particular energy and 
calculate Jic through equation Jic = — E). 

Then we choose a pair of {J, Jz), J can be a few times 
larger than Jjc but Jz keeps smaller than Ji^. Given the 
combination of {E,J,Jz) and the potential distribution 
we can find the apocenter position given by (r, 9) . Here r 
is distance to center and 9 is the angle between position 
vector and z-axis. We note that there are actually four 
parameters {E, J, Jz,9) to define the initial conditions 
for a particular orbit. So we further sampled 100 data 
points in 9 dimension. In order to plot the result in a 2D 
plane, we introduce a filling factor P for every (E, J, Jz) 
combination to describe this 9 dependence, which is the 
fraction of stars in the loss cone for a given combination 
of {E, J, Jz) (number of data points in loss cone divided 
by total sample size, e.g. 100 ), meaning that among all 
stars with same {E, J, Jz) only a fraction of P are inside 
the loss cone. 

By our definition a star in the loss cone will be dis¬ 
rupted by the BH within one dynamical time, so for ev¬ 
ery test particle we only integrate their orbits for one 
orbital cycle. If a particle comes back to its apocenter, 
we consider it as out of the loss cone and move to the 
next integration with new initial orbital data. 

Fig. |4] shows results from the experiment in a slowly 
rotating model (wq = 0.3), it represents the loss cone 
shape in phase space. Since J is not conserved we use 
its initial value at apocenter for the figure; at the time 
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Figure 4. Filling factor P of the loss coneas a function of J and Jz 
for different energies; the x axis is the modulus of the total angular 
momentum J, while the y axis is its z component Jz- Panel a) 
correspond to E = —1.3, b) = —1.5, c) E = —1.7, d) E = —1.9. 
Colors represent the filling factor P in percentage. All data are 
given in standard A-body units. 

of disruption J must be less than Jic. From panel a) 
to d), the energy of the test particles are in descend¬ 
ing sequence, so their position of apocenters are getting 
closer and closer to central BH. One can see that the 
whole plane comprises 3 regions: 1) inner region where 
P equals 1, meaning particles with these {E,J,Jz) can 
hit the BH within one dynamical time scale; 2) transi¬ 
tion region where P is non-zero but less than 1, particles 
with these {E, J, Jz) have a chance to hit the BH depend¬ 
ing on their apocenter position {9 value); 3) outer region 
where P = 0, none of particle in this region can hit the 
BH. In panel a) one can see only a few points are red and 
a lot of points are located in transition region. From a) 
to d), the fraction of P = 1 points in the (J, Jz)-plane 
increases and the transition region is compressed by the 
inner and outer region in horizontal direction (J dimen¬ 
sion). This is because test particles with high energy 
(loosely bound or unbound with respect to the BH) can 
go beyond the BH’s influence radius to the intermediate 
and outer regions of the cluster, where the axisymmetric 
stellar potential dominates. The angular momentum of 
these test particles will have large variations. So a wide 
transition region exists in high energy cases. But in the 
low energy case (stars strongly bound to the BH), e.g. 
panel d), test particles are moving inside the BH’s in¬ 
fluence sphere where the potential is dominated by the 
BH and thus approaches spherical symmetry. All loss 
cone stars following the classical loss cone approxima¬ 
tion, should have both J and Jz to be smaller than Jic. 
In all panels of Fig. IH on the contrary, we see how stars 
with J > Jic could be still in the new, extended loss 
cone of an axisymmetric system with a certain non-zero 
probability. 

For faster rotating models (wq = 0.6) the results are 
similar. Three regions are presented on the {J, Jz) plane, 
however, the extent of each region is different from the 
counterpart of same energy in slow rotating model. Fig. 
El gives an example, in both left (wq = 0.3) and right 
(ojQ = 0.6) panel the test particle have same energy, 
however, the resulting appearances are quite different. 
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Figure 5. Comparison between 2 models with same test particle 
energy —1.5. x axis is module of angular momentum in A^-body 
unit, y axis is z component of angular momentum. Left panel 
correspond to lu = 0.3; right panel uj = 0.6. Colors indicate the 
filling factor P in percentage. 

In the left panel we see the the outer border extended to 
J = 0.024, while in the right one the outer border goes 
to J = 0.04 and is not as clear as that in the left panel. 
Also in the right panel the red region is almost disap¬ 
peared. These results show how rotation modifies the 
loss cone shape in phase space. In both of these plots, 
the maximum radius stars can achieve are roughly the 
same. However, faster rotation means we have a more 
flattened cluster shape, which enhances the torque act¬ 
ing on stars, thus the variation in J becomes larger. So, 
the higher the degree of rotation in the stellar system, 
the larger is the extension of the loss cone in J direction. 

On first glance at Fig. |4] (also Fig. [5]) one might think 
that the loss cone is generally enlarged by a significant 
factor. However, as we pointed out above, there is a 
filling factor P for every point on the {J,Jz) plane. To 
find the net enlargement of the loss cone in axisymmetric 
potentials we introduce an effective area S of the loss 
cone in these plots by integrating the filling factor P over 
the (J, Jz) plane. For example, the effective area of the 
classical loss cone is just given by the size of the triangle 
J, Jz < Jic in our plots, since in the classical case P is 
unity everywhere in this triangle region. 

Now we compare the loss-cone size comparing the 
integrals S with each other. We define the quotient 
c^ic = Seff/Sic, where Sic = J’ld"^ is the classical loss 
cone integral, aic is plotted in Fig. [6] as a function of 
binding energy \E\. In the plot we show both slow and 
fast rotating models at two different evolution times. For 
the slow rotating model the ratio aic is even smaller than 
unity for binding energies larger than 1.4 - 1.5, meaning 
that at large \E\ the loss cone is smaller compared to 
the classical one. This is caused by the reduction of the 
probability P at the boundary and inside the classical 
loss cone region J < Jic- P is decreasing from inside to¬ 
ward outside. While for intermediate \E\ case, although 
P is still decreasing function of J, the large number of 
valid points overwhelms, so the net effect is increasing 
the effective area. However, if one goes further toward 
small \E\ the ratio will drop again, like the case of fast ro¬ 
tating model. This is just because P is sufficiently small 
in this case and win the game. For fast rotating model, 
another interesting feature is the ratio drops below 1 at 
I if I = 1.8. From this figure we see that the enlargement 
of loss cone, quantified by the ratio aic as a function of 
binding energy. Interestingly, the change of the effective 
loss cone size in every energy slice is less than 5-10%. 
These mild changes seem unable to raise TDR with the 
amount observed in the simulation, to address this it will 
be useful if we can estimate the TDR based on the effec¬ 


tive loss cone measurement and compare with simulation. 
However, the knowledge of how stars are distributed in 
energy and angular momentum is required. With the 
limited particle numbers of the model cluster, it is diffi¬ 
cult to get an accurate and reliable distribution function. 
Also in current work, we sample the energy space with 
large intervals {AE = 0.1), which may cause large errors 
in the estimated TDR. So we did not make the estima¬ 
tion. There are still plenty of works could be done with 
this topic. 
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Figure 6. Ratio between effective area Seff of loss-cone in ax¬ 
isymmetric system and Sic iri spherical system. 


5. ORBITAL PROPERTIES OF DISRUPTED STARS 

In this section, we investigate the origin of disrupted 
stars. Under the assumption that stars in loss cone can 
survive for only one orbital period, the origin of these 
stars can be examined by looking at their energy and 
angular momentum, as well as their origin (apocenter) 
in spatial coordinates (radius and angle 9). In spherical 
systems one can use effective potential to compute the 
apocenter of orbit, but in the axisymmetric case we do 
not have such convenient solutions except to run the sim¬ 
ulation twice. In the first run we find out the ID for those 
stars that will be disrupted by BH. Then, in the second 
run we make records for these stars more frequently than 
other stars, in order to catch their last apocenter posi¬ 
tion. 

We found in the beginning that the total TDR, espe¬ 
cially for small rt does only marginally depend on the 
rotation of the system; consistent with this we found in 
the previous chapter that the loss cone structure does 
change significantly, but the total integral over the loss 
cone space for axisymmetric systems yields only rela¬ 
tively small changes. Still it is interesting to study how 
the orbital properties of stars, which are tidally dis¬ 
rupted, change with the rotation of the system. In order 
to address this, we now turn back to our full N body 
simulations and study the distribution of \E\ (Fig. [7| 
and J (Fig. [8]) of the disrupted stars at their apocenter 
passage in three time intervals. From Fig. [7] one can see 
that most of the tidally disrupted stars have a binding 
energy between 1 and 2, coincident with the small bumps 
in Fig. [6] where aic > 1- Another evidence comes from 
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Figure 7. Panel a) and b) show normalized distribution of binding 
energy \E\ of tidally disrupted stars, for different rotating models 
and for three different time intervals (indicated by color) in the full 
A^-body simulation with rt = 10“^. The distribution is normalized 
to the total number of disrupted stars in each time interval. Panel 
c) and d) show cumulative fraction profile corresponding to a) and 
b), respectively. 
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Figure 8. Same as Fig. [T] but here for the distribution of total 
angular momentum of the disrupted stars. 


the distribution of J as shown in Fig. [8l where one can 
see the peaks are lying outside of the J/c which is roughly 
0.015. The peaks are moving toward larger J, which is 
caused by the increase of BH mass (recall the expression 
for Jic). A significant fraction of stars comes from places 
outside of the classical loss cone in (J, J^) plane. 

In spherical systems it is usually sufficient to describe 
the apocenter of an orbit by its radial distance from the 
center (the BH); the orientation of the orbit does not play 
any role for the orbital time and the nature of the en¬ 
counter with the central BH. However, in axisymmetric 
systems, orbits with diflFerent angle 9 (the angle between 
position vector of the star at apocenter and the z-axis) 
will differ from each other significantly. Therefore we 
have to describe the distribution of apocenters of tidally 
disrupted stars in terms of both the r (Fig. [5]) and 6 
(Fig. [TD]) dimension. From Fig. [5] one can see the peaks 
are quite far from the BH, in a region comparable to the 
BH influence radius, which is similar to the apocenter 



R [NB] 


Figure 9. Distribution of last apocenter distance of disrupted 
stars in 3 different time interval. Each curve are normalized to the 
total number of disrupted stars in given time interval. 




Figure 10. Normalized distribution of zenith angle 8 of last apoc¬ 
enter of disrupted stars. Left panel is axisymmetric model, right 
panel is spherical model. 

distribution in spherical systems (Paper I). The differ¬ 
ence turns out to be in the 9 dimension, as shown in 
Fig. [TUI We compare the 9 distribution between spher¬ 
ical and axisymmetric systems. Imagine we project all 
the apocenter points onto a sphere with radii equals 1. 
The measured number counts in each 9 bin AN{9) are 
computed by 27r • S(9) ■ sm(9)A9, where S(9) is the sur¬ 
face density of projected points on the unit sphere. If 
apocenters are uniformly distributed with 9, S(9) is con¬ 
stant, then AN(9) oc sm{9)A9. Here we choose an equal 
bin size, so the measured number count should follow a 
sin(0) curve. The right panel of Fig. |TU| plots 9 distri¬ 
bution for spherical model, which is taken from our last 
work (Paper I). In left panel we see the last apocenter dis¬ 
tribution have deficit at polar region comparing to sin(0) 
curve, and excess at places beyond and below the equato¬ 
rial plane, showing a double peak feature. The deficit at 
the polar region may have something to do with the flat¬ 
tening of the cluster, however, this is not the only reason. 
The double peak feature around the equatorial plane ob¬ 
viously does not relate to a geometrical origin, otherwise 
the peak should be placed at the equatorial plane. In 
Fig. 111! we compare the 9 distribution between slow and 
fast rotating models. One can see that in fast rotating 
model, the double peaks are more significant, accompa¬ 
nied by a further deficit in the angle range 0.2 — O.Ttt and 
0 .6-0.87r. 

In order to understand the double peak feature, we 
turn to the orbit structure of these disrupted stars. In 
non-spherical symmetric stellar system with a SMBH 
in its center, the space populated by stars can be di¬ 
vided into three parts depending on the distance to 
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Figure 11. Normalized distribution of zenith angle B of last apoc- 
enter of disrupted stars. Compare between slow (w030) and fast 
(w060) rotating models. 

the BH, namely the re gular, chaotic and mixing region 
(|Poon fc Merrittil2001f) . Inside the BH’s influence radius 
rfi, the potential felt by the star is dominated by the BH 
plus a small perturbation from the non-spherical stellar 
potential. In this region, the motion of stars is essen¬ 
tially regular, as in a spherical potential. Outside of rh, 
stars passing the center will suffer a large angle deflection 
by the BH, which in conjunction with the non-spherical 
potential near and outside Vh, could make their orbits 
stochastic. 

We are interested in stellar orbits in an axisymmet- 
ric stellar potential, which can get close to the central 
BH. These are typically two classes of orbits , short-axis 
tube (SAT) and saucer (see iVasilievI (j20I4f ) for exam¬ 
ple); they can be distinguished by their third integral 
of motion I 3 . Although I 3 may help us quickly distin¬ 
guish orbit f amilies, hnding the func tional form of I 3 is 
difficult fsee iLupton &: GunnI (|I987ll and discussion in 
iSridhar fc Toumal (j 1 999( 11 and is beyond the scope of this 
paper. We choose alternative ways to do orbit classihca- 
tion, such as Surface of Section (SoS) plot and Fourier 
analysis of Jx (see Appendix). 



due to restrictions by the 3rd integral. We also check the 
value of J at each apocenter passage. We find that SAT 
orbit achieve its minimum J at the equatorial plane; a 
saucer orbit cannot reach the equatorial plane, but its 
minimum J is achieved at the place which is next to the 
equatorial plane as marked in the plot by A plane. Re¬ 
call in the last section we said no matter what J one star 
has at the apocenter, at the time of disruption it must 
be smaller than Jic- So the last apocenter place should 
be around the A plane. This seems to be promising to 
explain the double peak in 9 distribution, however, need 
to be confirmed. In order to see this we try to do orbit 
classification for the disrupted stars, which is computa¬ 
tionally expensive. So we just randomly select a sub¬ 
sample of disrupted stars and divide them into 3 orbit 
families: SAT, saucer and others (here “others” means 
they do not belong to the former two families, and may 
be chaotic orbits). Among the 2943 sample stars, 1719 
are classified as “others”, 757 as saucer and 467 as SAT. 
Then we re-plot the r and 9 distribution for different 
orbit families in Fig. [131 
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Figure 12. Examples of orbital structure for SAT (left) and 
saucer (right) orbit. Stars achieve maximum J when their instant 
orbital plane coincide with B plane, while minimum when coincide 
with A plane. 

Fig. |T2] gives examples of SAT and saucer orbits in 
configuration space. The plot is made in cylindrical co¬ 
ordinates so that one can catch the main point easily. 
For SAT orbit, one can see its apocenter can go both 
above and below the equatorial plane. While apocenter 
of saucer orbit can only exist on one side of mid-plane, 


Figure 13. Normalized distribution of apocenter distance r and 
zenith angle 9 of last apocenter of disrupted stars for different orbit 
families. The distribution is normalized to the number of stars in 
each orbit family. 

The results show that the apocenter distribution of dif¬ 
ferent orbit families not only differs in 9 but also in r. 
One can see the innermost region is dominated by SAT 
orbits, and concentrated to the equatorial plane. Inter¬ 
mediate radius is mostly occupied by saucer orbits, and 
the distribution in 9 shows double peaks as expected. 
Further out is the region dominated by orbits marked 
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as others. These orbits can go outside of the influence 
radius and are basically chaotic orbits. From Fig. [13] 
one can also find out the fractions of each orbit family 
contributing to the budget of disrupted stars: the largest 
fraction comes from chaotic orbits; SAT orbits contribute 
least to the budget because they are deeply buried in 
the cluster center where the total star number is small; 
the intermediate contribution is from saucer orbits which 
create the two peaks in the 9 distribution. 

6. CONCLUSIONS AND DISCUSSIONS 

Tidal Disruption (TD) of stars by supermassive central 
black holes (SMBH) from dense rotating star clusters is 
modelled by high-accuracy direct iV-body simulation. As 
in a previous paper on spherical star clusters we study the 
time evolution of the stellar tidal disruption rate and the 
origin of tidally disrupted stars, now according to several 
classes of orbits which only occur in rotating axisymmet- 
ric systems (short axis tube and saucer). In empty loss 
cone regime, comparing spherically symmetric and ax- 
isymmetric systems we find a higher TD rate in large 
rt models in axisymmetric case, but for small rt case - 
somewhat surprisingly - there is virtually no difference in 
the TD rate, maybe a small increase due to axisymmetry. 

We define an extended loss cone by the condition that 
stars in the axisymmetric potential reach the BH within 
one orbit. A detailed analysis shows that the structure 
of the loss cone significantly differs from the spherical 
case; if Jic is the critical angular momentum to be in the 
loss cone in a spherical system, and J, Jz are the total 
and z-component of the angular momentum of a stellar 
orbit, there are many stars with J > Jic in the loss cone; 
since, however, there are also some stars with J > Jic, 
which are now not in the loss cone. In the total balance 
the number of loss cone stars is only very moderately 
increased. 

In the experiment of measuring the shape of loss cone, 
we assume the test star can survive only one dynami¬ 
cal time in collisional system, after one orbit it will be 
“kicked” to another place in phase space due to interac¬ 
tions with other stars. However, in collisionless limit, if 
we allow the test star to survive more orbit cycles, test 
star with much higher J will also have chance to get rid of 
its angular momentum and be disrupted by BH. Then it 
is possible that an even larger loss cone region in phase 
space than what we presented here may exist, and re¬ 
sult in a higher TDR. In order to check this, simulations 
with much more particles are needed and we would like 
to leave this task for future work. 

The orbit type of disrupted stars strongly depends on 
energy as we discuss in detail in the previous sections. 
TD of stars most strongly bound to the BH are domi¬ 
nated by short-axis tube (SAT) orbits. In intermediate 
regions saucer orbits dominate, which create a character¬ 
istic double peak structure in the last apocenter position 
of their orbit relative to the equatorial plane. And fur¬ 
ther out chaotic orbits. 

It is known for almost half a century that tidal 
disruption of stars should occur near SMBH, but 
only much more recently the X-ray e mission o f tidal 
disruption events has be en detected ([Komossal 120021 : 
iKomossa &: Merritt! 1200^ . A sim ple a. r gumen t on the 
fallback time for tidal debris by iReesI (1198811 has led 
to the prediction of a characteristic power law of the 


light curve with time, which can be used to distinguish 
TD events from other transients. It is interesting that 
a SMBH binary can cause characteristic disr uptions in 
such an otherwise standard TD light curve (ILiu et al.l 
1201411 . Hayasaki and collaborators claim that eccentric 
TD events lead to somewhat longer liv ed central accre¬ 
tion disks ([Havasaki et al.l I2013L1201^ . It will be very 
interesting to see whether and how the evolution of tidal 
debris and the fallback rate are affected by different or¬ 
bits of the disrupted stars as discussed here. 

It has been claimed that rotation may help to quickly 
refill loss cones around binary supermassive black holes, 
which helps significantly to accelerate shrinking and 
final coalescence of SMBH binaries in cosmologically 
short time scales (iBerczik et al.ll2005 iPreto et al.l 20111 : 
IKhan et'al]l2013t IKhanll2014ll . In our paper we study by 
direct A^-body simulation the tidal accretion of stars and 
their orbital parameters in r otating axisymmetr i c sys - 
tems. We confirm the result of iVasiliev &: MerrittI (|2013ll 
that there is an increase in the rate of refilling of the loss 
cone, but it is moderate. However, the situation deserves 
more detailed study, because a SMBH binary creates a 
much stronger deviation from spherical symmetry than 
the one used in our models with single SMBH. And sec¬ 
ond the detailed structure of the rotation in a central 
nuclear star cluster could affect the enhancement of the 
loss cone. 
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APPENDIX 

ORBIT CLASSIFICATION 
Surface of Section 

From Fig.[T3]we can see the whole accessible region on (i?, vr) plane is divided into two parts (note that points with 
opposite VR actually belongs to same orbit, so this plot is symmetric with horizontal axis). Each part represents a 
family of orbit. Curves that intersect with R-axis are footprints of short axis tube (SAT) orbits, others are of saucer 
orbits. 



Figure 14. Surface of section plot, x axis is distance R to origin on equatorial plane, y axis is vr = dR/dt when star go across the 
equatorial plane. 


Fourier Analysis of evolution 

In axisymmetric potential, force is not centripetal hence exerted a torque on the star which will change the x and 
y components of its angular momentum. Fig 1151 show the time evolution of Jx for both SAT and saucer orbits. The 
pattern of Jx and Jy are the same but shifted with a phase of 7r/2, so in the following discussion we only focus on Jx- 
Furthermore, the evolution of Jx shows some quasi-periodicity. From eye inspection, one can guess the mathematical 
expressions for the curves. 

As shown in Fig. [T5l the curve for SAT orbits seems to be represented by sin(fit)(l -I- cos(f 2 t)) (fi < f 2 ), which 
can be further converted to sin(fat) -I- sin(fbt) -I- sin(fct) (ignore coefhcients before the trigonometric functions), with 
fa = fi, fb = f 2 — ff and fc = f 2 + fi ■ If we perform a Fourier analysis on this curve, we expect to find 3 principal 
frequencies: fa,fb,fc in ascending sequence. And these 3 frequencies satisfy the equation fa = (fc — fb)/2. 

For saucer orbits, the curve seems to be represented by sin(fit)cos(f 2 t) (fi < f 2 ), following the same procedure we 
expected to find 2 principal frequencies: fa, fb, with fa = f 2 — fi and fb = f 2 + fi- 

A demonstration is shown in Fig. 1161 one can clearly see the 3 principal frequencies for SAT orbits and the 2 for 
saucer orbits. Some of the small peaks appeared at higher frequencies which is the order harmonics and some are 
produced from other components. 

We use both methods to cross check the validity of orbit classification for the tidally disrupted stars. 
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Figure 15. Time evolution of Jx for SAT orbit (left) and saucer orbit (right). 
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Figure 16. Fourier frequency distribution of Jx- Horizontal axis is frequency in unit of [T] Vertical axis is amplitude of corresponding 
component. Left panel represents SAT orbit, right panel represents saucer orbit. 
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